arXiv:1503.07044v2 [quant-ph] 19 May 2015 


Multistable particle-field dynamics in cavity-generated optical lattices 


Dominik J. Winterauer/ Wolfgang Niedenzu,^’^ and Helmut Kitsch^’* 

^ Institut fiir Theoretische Physik, Universitdt Innsbruck, Technikerstrafie 25, A-6020 Innsbruck, Austria 
^Department of Chemieal Physies, Weizmann Institute of Science, Rehovot 7610001, Israel 

(Dated: May 19, 2015) 

Polarizable particles trapped in a resonator-sustained optical-lattice potential generate strong 
position-dependent backaction on the intracavity field. In the quantum regime, particles in different 
energy bands are connected to different intracavity light intensities and optical-lattice depths. This 
generates a highly nonlinear coupled particle-field dynamics. For a given pump strength and detun¬ 
ing, a factorizing mean-field approach predicts several self-consistent stationary solutions of strongly 
distinct photon numbers and motional states. Quantum Monte Carlo wave-function simulations of 
the master equation confirm these predictions and reveal complex multi-modal photon-number and 
particle-momentum distributions. Using larger nanoparticles in such a setup thus constitutes a 
well-controllable playground to study nonlinear quantum dynamics and the buildup of macroscopic 
quantum superpositions at the quantum-classical boundary. 

PACS numbers: 37.30.+i,37.10.Vz 


I. INTRODUCTION 

The dynamics of polarizable point-like particles 
trapped by an optical cavity light field has been the sub¬ 
ject of intense theoretical and experimental studies in 
the past decade [1-4]. Beyond implementing improved 
neutral-atom cavity-QED systems [5-7], recently pro¬ 
posed applications of such setups range from ultrahigh- 
Q optomechanics [3, 8, 9] to precision tests of quantum 
mechanics at a mesoscopic scale [10] and gravity [11[. 
Following the first pioneering experiments more than a 
decade ago [5-7], several groups have implemented reli¬ 
able cavity-based optical traps in their experiments for 
various particle numbers ranging from a single or few 
atoms [12-15] to Bose-Einstein condensates (BECs) [16- 
19] or lately even considerably heavier nanoparticles [20- 
22 ]. 

Particles in cavity fields, in contrast to free-space opti¬ 
cal potentials, substantially act back on the field dynam¬ 
ics [23, 24], which generates complex and rich nonlinear 
dynamics [25-28]. In the standard optomechanical limit 
of very tightly trapped particles or membranes, which can 
essentially be modeled by harmonic oscillators [29, 30], a 
wealth of interesting physics beyond ground-state cooling 
appears in the strong-coupling regime. Typical examples 
are atom-field entanglement, nonlinear oscillations, and 
multistable behavior [26, 31-35]. The system dynamics 
gets even more complex and rich, if one refrains from 
linearizing the particle motion and considers its full dy¬ 
namics along the cavity potential [36, 37] . 

In most cases the optical potential along the cavity axis 
is well approximated by a sinusoidal lattice potential with 
a depth proportional to the momentary intracavity pho¬ 
ton number [2[. While for deep potentials the harmonic- 
oscillator basis allows for analytic insight, it becomes in¬ 
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adequate for shallower lattices. The eigenfunctions of pe¬ 
riodic potentials are delocalized Bloch functions, which 
can be transformed to localized Wannier functions [38] . 
Unfortunately, no analytic solutions for neither the Bloch 
nor the Wannier functions are known even for a fixed lat¬ 
tice depth. Hence, aiming for an explicit analytic treat¬ 
ment in the (dynamic) quantum-potential limit is a hope¬ 
less goal. In view of these complications, several semiclas- 
sical and mean-field models with factorized evolution of 
the particles and the field have been developed to obtain 
some first insights [23, 26, 39] . Here the field expecta¬ 
tion value is governed by ordinary differential equations 
containing particle expectation values. This field is in 
turn inserted in the effective Hamiltonian for the parti¬ 
cle motion [36, 40] . Even in this strongly simplified limit 
the nonlinearity of the interaction does not allow for a 
straightforward solution in the general case and further 
assumptions are needed [33] . 

In this paper we study the full quantum dynamics and 
the steady-state properties for the case of a single par¬ 
ticle in a cavity-sustained optical lattice in the strongly 
coupled and strongly pumped limit. Hence, our treat¬ 
ment will centrally be based on straightforward numeri¬ 
cal solutions of the corresponding quantum-optical mas¬ 
ter equation. Strong emphasis will be put on steady-state 
properties of the system in the limit of very low tempera¬ 
tures close to T = 0, where semiclassical treatments pre¬ 
dict a multitude of stationary solutions. To this end we 
will heavily rely on quantum Monte Carlo wave-function 
simulations [41-43], since a direct solution of the master 
equation becomes very slow and cumbersome owing to 
the large joint particle-field Hilbert space, even though 
we consider the simplest possible system involving only 
a single particle. 

This paper is organized as follows. In Sec. H we intro¬ 
duce the model Hamiltonian and the master equation, 
from which we derive equations for the expectation val¬ 
ues of the cavity field and the photon number, depending 
on the particle state. To get some first qualitative in- 
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Figure 1. (Color online) A particle within a driven optical 
cavity. The longitudinal cavity pump rj builds up an intra¬ 
cavity field that drives the particle motion. The particle’s 
motional state affects the cavity detuning (dynamic refrac¬ 
tive index), which in turn influences the intracavity photon 
number. Photons leak through the cavity mirrors at a rate 
2k. 

sight into the system behavior and to identify interesting 
parameter regimes, we start with simplified semiclassical 
models. Factorizing atom and field dynamics, we approx¬ 
imate the photon field by a classical field characterized 
by its mean photon number, which is determined by the 
spatial distribution of the particle. We then look for self- 
consistent steady-state solutions for the expected photon 
number. Section III is devoted to studies of these self- 
consistency conditions in various limiting cases. We first 
consider the deep-trap limit of harmonic particle confine¬ 
ment, which allows for an analytic treatment. This anal¬ 
ysis is afterwards extended to localized Wannier states in 
the general case of a periodic optical lattice. In Sec. IV we 
then numerically solve the full master equation in typi¬ 
cal operating regimes determined before and also analyze 
the behavior of single Monte Carlo trajectories. Finally, 
in Sec. V the conclusions are drawn. 


II. MODEL 

We consider the standard case of a driven, damped 
cavity mode with a single polarizable particle of mass 
m moving along the cavity axis, as sketched in Fig. 1. 
The Hamiltonian (fi. = 1) in a rotating frame with pump 
amplitude t], cavity detuning Aq, and effective particle- 
field interaction strength Uq is then given by [44] 

i? = ^ - [Ac - Uo cos^(A:Ra:)] a^a -ir][a- a'^') , (1) 

where = 27r/A is the single-photon recoil momentum, 
with A being the cavity mode wavelength. The particle 
position and momentum operators x and p, and the pho¬ 
ton mode annihilation (creation) operators obey the 
standard canonical commutation relations 

[x, p] = i, [a, a'l'] = 1, (2) 

with all other commutators vanishing. The coupling 
strength is parametrized by Uq, which denotes the op¬ 
tical potential depth per photon as well as the maximum 
cavity mode resonance frequency shift that a particle in¬ 
duces when placed at a field antinode. Here we consider a 


large negative Uq, which corresponds to high-field seeking 
particles. 

We assume operation far from any internal optical res¬ 
onances, such that spontaneous light scattering and ab¬ 
sorption losses from the particle into the mode can be 
neglected [1]. The dominant loss mechanism is then 
cavity damping, which at optical frequencies can be in¬ 
corporated by a standard master equation treatment 
parametrized by a loss rate k [45] , 

p = —i [H, p] + K (2apa^ — a^ap — pa'a) . (3) 

From this master equation we straightforwardly derive 
ordinary differential equations (ODEs) for the expecta¬ 
tion values of the field amplitude, (a) = a, and the pho¬ 
ton number, (n) = (a^a), which read 

a = [i (Ac — Uq (cos^(fcRa;))) — k] a + p (4a) 
(n) = p{a + a*) — 2 k (n) . (4b) 

Within the semiclassical treatment with a c-number de¬ 
scription of the field amplitude the particle-field density 
matrix is assumed to be separable. Obviously, the field 
dynamics depends on the motional state of the particle 
via the expectation value (bunching parameter) 

b := (cos^{k^^x)) (5) 

in a nonlinear fashion. This parameter itself is, in 
turn, governed by the Schrodinger equation contain¬ 
ing the Hamiltonian (1), whose spatial eigenstates are 
Bloch functions according to the quantized lattice depth 
Vq = |f7o|a^a. This yields a different evolution for each 
photon-number component of the total wave function and 
thus a very complex time evolution. Hence, a full solu¬ 
tion of the master equation (3) requires a numerical ap¬ 
proach, which can be directly implemented using a trun¬ 
cated photon number and momentum basis expansion. 
Note that due to the periodic nature of the potential we 
can work with periodic boundary conditions in real space 
and use a discrete momentum basis {jp) = |jT:r)} with 

j e No. 

As these calculations are time consuming and the range 
of physical parameters (p, Ac, Nq, wr) is large, we first 
try to get some qualitative insight and find interesting 
parameters regions using the factorized semiclassical ap¬ 
proach involving Eqs. (4). 

III. SELF-CONSISTENT SEMICLASSICAL 
SOLUTIONS OF THE COUPLED ATOM-FIELD 
DYNAMICS 

Let us now analyze potential stationary solutions of 
the coupled ODE system (4). As the field dynamics in 
the semiclassical approximation depends on the position 
distribution of the particle via the expectation value (5) 
only, for the system to reach a steady state we need a sta¬ 
tionary wave function. This leads to the self-consistency 
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condition 


+ (Ac - C/q (cos2(fcRa;)))^ ’ 

where the wave function of the particle has to be an 
eigenstate of Eq. (1) with the photon-number operator 
a^o replaced by (n). Note that the expectation value 
(cos^(fcRa;)) in the denominator on the right-hand side 
of Eq. (6) does not explicitly involve any field operators. 
Nevertheless, the time evolution of the spatial part of 
the wave function depends on the field intensity. Hence, 
the state can only be stationary if it is an eigenstate of 
the Hamiltonian (1) for the momentary photon number. 
Note that the pump amplitude ry is a free parameter in 
the above equation and in many cases for a given eigen¬ 
state of the particle Hamiltonian a self-consistent choice 
of rj can be made to fulfill Eq. (6) [40]. We, however, opt 
for the opposite approach and determine self-consistent 
photon numbers for given pump strengths. 

Let us mention, though, that this is only a necessary 
condition and by no means sufficient for a stable station¬ 
ary equilibrium subject to the quantum fluctuations of 
the system. At this point it can only serve as a guide 
towards interesting parameter regions, which is, e.g., the 
case when several different spatial eigenfunctions lead to 
the same pump amplitude rj. We will discuss this in some 
more detail below for specific limiting examples. 


A. Harmonic-oscillator expansion in a deep lattice 

In the limit where the potential depth Vq « |17o| (n) 
strongly exceeds the recoil energy Er = wr := k‘^/{2m), 
the lowest-energy particle states are well localized within 
a single well of the optical lattice. For low enough tem¬ 
peratures the optical potential Vee{x) = Uq {n) cos^(fcRx) 
can then be approximated by a harmonic potential [40], 

Uo (n) cos^(A:rx) « Uq (n) (l - . (7) 

The corresponding trapping frequency Who then reads 



Ac/wr 
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Figure 2. (Color online) Self-consistent photon-number con¬ 
tours for a harmonically trapped particle as a function of 
the cavity-pump detuning for different pump amplitudes rj. 
The four plots show contours in the Ac-{n) plane for dif¬ 
ferent harmonic-oscillator eigenstates [riho), where Eq. (11) 
holds self-consistently. All excited states, riho > 1, yield the 
possibility of two self-consistent solutions for a certain range 
of detuning Ac and therefore may allow for optomechanical 
bistability. Parameters: Uo = — IOwr and k = wr. 


Hence, within the harmonic-oscillator approximation 
Eq. (6) becomes the simple algebraic equation 


(n) = 


Ac - C/o 1 - 


2nho-H 


( 11 ) 


which can be easily solved for each choice of eigenstate 
number Uho- Figure 2 shows contours in the Ac-{n) plane 
for different values of rj for which Eq. (11) holds. While 
the lowest-energy state riho = 0 results in a unique pho¬ 
ton number, multiple (up to two) solutions are possible 
for higher excited states. 


^ho 



(n) > 1, 


B. Self-consistent states for the fnll lattice 
(8) dynamics 


and we can analytically find the respective oscillator 
states jriho) to this frequency. The expectation value in 
the denominator of Eq. (6) is then well approximated by 
(cos(/cRa;)^) « 1 — /cr with 

such that 

4{xb„.,. = (2"l.o + l)^ = ^=|fT, (10) 


Several interesting aspects of the deep-trap harmonic- 
oscillator regime have been studied in the past [27, 33, 
40] . In many respects the system is directly related 
to standard optomechanical models with quadratic cou¬ 
pling [29] . For ultracold particles and weaker optical 
potentials the motion is strongly delocalized in the lat¬ 
tice [2[. Hence, we now turn to the full model and con¬ 
sider the particle’s motion in the periodic optical lattice 

Ves{x) = Uo (n) cos^ik^ix). (12) 

For very shallow lattices close to zero temperature, i.e., 
for a BEC in a cavity, a two-mode expansion of the wave 
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function can be applied [16, 46], which again allows for 
analytic treatments and analogies with optomechanical 
couplings. However, the validity range of this model is 
limited in temperature, time, and coupling strength. As 
we are here more interested in the limit of strong nonlin¬ 
ear backaction in deep potentials, we cannot apply this 
simplification and have to solve the Schrodinger equation 
for a periodic potential, which gives us the well-known 
Bloch states '^rnq{x), where m denotes the energy band 
and q is the quasi-momentum [38]. Being periodic with 
the lattice constant, they are not the best basis to de¬ 
scribe a single localized particle. Hence, we switch to a 
Wannier basis, where each basis state represents a local¬ 
ized wave function with its center of mass at a particular 
lattice site. Such basis states have been very successfully 
used to study ultracold particle dynamics in optical lat¬ 
tices [47, 48] . The Wannier functions for a given band 
index m localized at lattice position R are defined as [38] 

Wm{x - R) ■■= J ^ (13) 


where a is the lattice periodicity. The Bloch functions 
'i’mqix) are only defined up to a phase. In order to ob¬ 
tain the maximally localized (i.e., real and exponentially 
decaying) Wannier functions, these phases need to be 
properly adjusted [38] . In what follows we choose for 
simplicity R = 0. 

We are now able to restate the self-consistency equa¬ 
tion ( 6 ) for each band index m as 


0 = 


- (n) , (14) 


with 


+ (Ac - Uob^y 

/ OO 

[u;m(a;)]^ cos^(fcRa;) dx. (15) 

-OO 


Contrary to the harmonic oscillator wave functions, there 
is no analytic expression for Wannier functions and we 
have to numerically solve the Schrodinger equation for 
each particular (n). Therefore (n) does not explicitly 
appear on the right-hand side of Eq. (14), but enters 
implicitly through the shape of the wave function. As 
before, we can obtain the contours where Eq. (14) holds 
self-consistently in the Ac-{n) plane for the same values 
of ? 7 ; see Fig. 3. 

The behaviors of the photon numbers for the lowest- 
energy states nho = 0 and m = 0 in both cases are 
very much alike, because the corresponding lowest bound 
states in both models are similar. Indeed, the maximally 
localized Wannier functions converge to the harmonic os¬ 
cillator functions for deep potentials [38] . For the higher- 
energy eigenstates, however, the photon numbers differ 
significantly. The reason behind is that in a harmonic 
oscillator all states are bound, while Wannier states for 
increasing m > 1 undergo a transition from bound to 
free states for a given photon number (i.e., potential 
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Figure 3. (Color online) Same as Fig. 2 for particles in lo¬ 
calized Wannier states. Photon numbers above the dashed 
lines for Wm>i correspond to bound Wannier states {E < 0). 
Higher bands exhibit up to three solutions of Eq. (14) for a 
given value of Ac. Parameters: Uq = — IOwr and k = 0 ;^. 


depth). Dashed lines in Fig. 3 indicate this boundary. 
For m > 1 sharp bends appear, yielding self-consistent 
contours reminiscent of nonlinear response curves. The 
origin of these peculiarity at the transition from free to 
bound states becomes evident, if we look at the spatial 
particle density of the respective Wannier states. Fig¬ 
ure 4 illustrates the behavior of the fourth band Wannier 
state W 4 for different mean photon numbers (i.e., poten¬ 
tial depths). The key quantity here is the expectation 
value of the bunching parameter [Eq. (15)[, which 
determines the backaction of the particle on the cavity 
field, i.e., its effective refractive index. For free parti¬ 
cles, (E) > 0, the wave function is barely localized and 
bm ^ y Around (E) « 0 the Wannier states localize 
around potential maxima, i.e., optical field nodes, which 
minimizes the backaction of the particle on the cavity, 
bm < y while for deeper potentials (E) < 0 and parti¬ 
cles are drawn towards field antinodes and the index of 
refraction increases with potential depth, bm> Thus 
the nonlinear behavior of the refractive index allows for 
multiple self-consistent solutions for certain ranges of the 
cavity detuning Ac- In particular, we also find solutions 
corresponding to unbound particle states (e.g., for 1^4 in 
Fig. 3). 


C. Stability of the self-consistent factorized 
solutions 

As we saw above, for certain parameter ranges in both 
the harmonic oscillator and the Wannier contour plots 
more than one self-consistent solution appears. Whether 
or not these solutions have significant relevance for the 
full system dynamics depends on their stability proper- 
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Figure 4. (Color online) Spatial particle distribution in a 
fourth band Wannier state for different photon numbers: (a) 
Free particle: The wave function is barely localized, 64 < 0.5. 
(b) Transition to a bound state: The wave function is localized 
at potential maxima, bi < 0.5. (c) Tight-binding regime: The 
wave function is strongly localized in a single potential well, 
hi > 0.5. 


ties, i.e., their response against small deviations in the 
photon number or the spatial distribution. Some quali¬ 
tative insight can already be gained by virtue of Eq. (6) . 
The right-hand side of Eq. (6) depends on the shape of 
the wave function in real space, which in our semiclassical 
model implicitly depends on (n). The term on the right- 
hand side of Eq. (6) determines the mean photon number 
that is allowed by the spatial part of the wave function in 
steady state. If it increases or decreases with (n) faster 
than (n) at a self-consistent point, one may assume that 
the self-consistent configuration is unstable. Therefore 
we find that at stable self-consistent configurations the 
inequality 
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Figure 5. (Color online) Contour plot for the right-hand side 
of Eq. (14) for the fourth band Wannier state. Solid (dotted) 
lines mark stable (unstable) self-consistent solutions accord¬ 
ing to Eq. (16). Parameters: r] = 6cjr, Uq = — IOojr, and 

ft = CJR. 


IV. NUMERICAL ANALYSIS OF THE FULL 
COUPLED ATOM-FIELD DYNAMICS 


In order to test the above analysis, we now strive to 
solve the full master equation [Eq. (3)]. As already men¬ 
tioned, the (even for a single particle) large Hilbert space 
makes a direct numerical integration attempt practically 
unfeasible for realistic parameters and photon numbers. 
We therefore make use of quantum Monte Carlo wave- 
function simulations, in which single stochastic state vec¬ 
tors (instead of the whole density matrix) are evolved 
subject to a non-Hermitian effective Hamilton opera¬ 
tor [41-43]. This evolution is stochastically interrupted 
by quantum jumps corresponding to a projective removal 
of one photon. Mathematically, averaging over a large 
number of such trajectories then approximates the full 
density matrix. Interestingly, the jumps can be physi¬ 
cally interpreted as detection events of photons leaking 
out of the resonator. Hence, single trajectories provide 
a microscopic view of the processes incorporated in the 
master equation since the ensemble average over many 
trajectories converges towards the solution of the latter. 


d rf 

d (n) /c 2 _|_ (cos2(fcRa;)))^ 


- 1 < 0 


(16) 


must hold. This rather intuitive result is verified in 
Appendix A via linear stability analysis. The stabil¬ 
ity regions for the fourth band (where up to three self- 
consistent solutions exist) are shown in Fig. 5. 


In what follows we compare our above predictions 
with the time evolution of single trajectories as well as 
to their ensemble average. The numerical implementa¬ 
tion of these simulations was done within the C++QED 
framework allowing efficient and fast simulations [49-51] . 
Since dynamic aspects have been eliminated in our self- 
consistency and stability analysis, it is not clear which of 
the self-consistent solutions appear in the dynamics and 
what are their corresponding probabilities. 
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A. Time evolution of single trajectories 

We consider a small sample of single quantnm trajec¬ 
tories in a multistable regime. Figure 6 shows the cor¬ 
responding expectation values of the intracavity photon 
number (n) as well as of the kinetic energy (i?kin) = 
(p2) /(2m). As one might expect, both quantities jump 
simultaneously between rather stable values. The lat¬ 
ter can be identified as the possible semiclassical values 
found above. Each trajectory thus seems to switch be¬ 
tween these states rather than forming state superposi¬ 
tions. Between jumps both quantities appear to fluctuate 
only weakly about the self-consistent values (upper three 
graphs). In some cases (n) jumps to very low values, 
where no bound state exists. In such cases the system 
continuously heats up (i.e., (Akin) increases) until a sub¬ 
sequent jump occurs and projects the particle back into 
a bound state (as for example in Fig. 6b between the 
two quantum jumps at wnt « 120 and 150; the signif¬ 
icantly increased photon number after the second jump 
allows again for bound states). Figure 6d shows an ex¬ 
treme case, where the particle remains essentially free for 
a long time. We find that there exists a multitude of sta¬ 
ble self-consistent solutions of Eq. (14) around (n) = 4 for 
higher bands (m > 6), whose self-consistent photon num¬ 
bers increase only slightly with increasing band index. A 
small plateau of (Akin) in Fig. 6d can be interpreted as 
an occupation of the 12th Wannier state, (Akin)m=i 2 ' 
We also indicate the mean kinetic energy of the 14th ex¬ 
cited band, (Akin),n^i 4 , and observe that this energy is 
reached continuously rather than by discrete jumps as in 
the bound case. Below we will re-encounter reminiscences 
of such trajectories in the ensemble-averaged solution of 
the master equation. 

Due to the self-consistent photon numbers’ small sen¬ 
sitivity on the band index for m > 6, according to our 
semiclassical analysis several excited bands can co-exist 
at a given value of (n). Wannier functions for a given po¬ 
tential depth are mutually orthogonal, therefore transi¬ 
tions between bands can only occur through fluctuations 
in the unbound regime, yielding much slower transition 
rates than in the bound regime. Trajectories may jump 
back to bound states. Figure 6a-c, or remain unbound. 
Figure 6d. The likeliness of jumping back to a bound 
state seems to decrease with kinetic energy. At this stage 
it appears that the momentum part of the wave function 
controls the expected intracavity photon number rather 
than vice versa. 

The correlated particle-field jumps reflect strong 
particle-field correlations and some amount of entangle¬ 
ment, as previously discussed in similar contexts [34, 52]. 


B. Stationary solution of the master equation via 
ensemble-averaged quantum trajectories 

We now investigate the solution of the master equa¬ 
tion (3) of the joint particle-field density matrix by av- 
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Figure 6. (Color online) Expectation values on single quan¬ 
tum trajectories for rj = 6cjr, Ac = — T.Swr, Uq = — IOcur, 
and K = wr. The mean values jump between the stable val¬ 
ues predicted by the self-consistency equation (14) (dashed 
lines), (a) and (b) Trajectories with several jumps from very 
high to very low photon numbers, (c) Trajectory with only 
a few jumps, (d) Trajectory with a large increase of (Akin) 
during evolution in a low-photon-number state. The kinetic 
energy mean value of the 12th and 14th excited bands are 
shown as black solid lines. Although numerical values indi¬ 
cate that the trajectory could be in a definite band, the neat 
picture of correlated photon number and momentum jumps 
clearly breaks down at this point. Note the different scaling 
of the (Akin) axis. 


eraging over a sufficiently large ensemble of Monte Carlo 
trajectories. First we check the distribution of photon 
numbers for a specific choice of parameters and compare 
it to the semiclassical results. In Fig. 7 we depict the sim¬ 
plest case of parameters, where only a single semiclassical 
solution exists. Interestingly, we see that the mean pho¬ 
ton numbers obtained from the Monte Carlo simulations 
agree surprisingly well with the self-consistent solutions 
of Eq. (14) for the lowest Wannier state, as long as the 
cooling regime (large negative effective detuning) is main¬ 
tained. Closer to resonance we see a deviation towards 
higher photon numbers, which indicates the appearance 
of motional exited states (cf. Fig. 3). 
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Figure 7. (Color online) Self-consistent solutions of Eq. (14) 
for the lowest Wannier states (m = 0) compared to ensemble- 
averaged quantum simulation results. Solid lines mark self- 
consistent contours; data points are quantum simulation re¬ 
sults. The black line separates the regions where Aeff,o < 
— Aeff ,2 (left) and Aeff.o > — Aeg, 2 (right) along the self- 
consistent contours of the zeroth band; cf. Eq. (18). In the 
left region quantum simulations are in accordance with the 
self-consistent solutions of Eq. (14) for the lowest band, while 
deviations arise in the right region due to population of the 
m = 2 Wannier state. Parameters: Uo = —IOcur and k = cur. 

Motivated by the single trajectories depicted in Fig. 6 
one can deduce a microscopic interpretation of the dy¬ 
namics. To each band may be assigned a kinetic temper¬ 
ature ksT = 2Ekin = (p^) /to, which increases with the 
band index [53]. The sign of the effective detuning 

Aeff,m := Ac — C/o^m (17) 

determines whether the according band is heated (+) or 
cooled (—). Heating means that in a certain band the sys¬ 
tem tends towards populating higher excited bands, while 
cooling implies the opposite. Since the value of Aeff, m is 
different for every band, some bands (the lower ones) are 
heated and others (the higher ones) are cooled. From the 
proportionality of the cooling/heating rates to Aes, m, 
we conclude that higher bands appear in the ensemble- 
averaged steady-state solution if 

Aeff, m > — Aeff, m+2- (18) 

Note that for symmetry reasons the dynamics induced 
by the Hamiltonian (1) conserves the parity of the initial 
state. For a particle initially in the ground state, the low¬ 
est accessible excited state is the second band and conse¬ 
quently to-I- 2 appears in Eq. (18). Hence the system effec¬ 
tively remains in the lowest band until Aes^o > ~Aoff, 2 ; 
see Figs. 7 and 8. This implies that, though the system 
is effectively blue detuned, it does not get heated; see 
Fig. 9. For certain parameter values a further increase of 
Ac around Aes, o = 0 even yields further cooling before 
the second excited band is populated. 

A more complete picture of this dynamics can be found 
if one includes momentum distributions, as depicted in 



Figure 8. (Color online) Combined photon-number and 
momentum-state occupation probability after a time interval 
of At = starting from a state with zero momentum 

and one cavity photon, fc = 0 and n = 1, for Uo = — IOcur, 
K = u;r, rj = 6cjr, and different values of the detuning. The 
density matrix is approximated via quantum simulations with 
ensemble averages over 1000 trajectories for each parameter 
set. The numbers in each box give the detuning Ac in units 
of u;r. 



Figure 9. (Color online) Quantum simulation results show¬ 
ing the particle’s average kinetic energy (Akin) as a measure 
of its temperature as a function of the cavity detuning for 
different pump amplitudes. For rj = Scur an interval where 
the hnal temperature decreases with increasing detuning is 
clearly visible. The plot shows only a detail of the set of data 
points for better visibility of the effect. In order to elimi¬ 
nate temporal fluctuations at single time instances the plot¬ 
ted points are time averages of (Akin) of 100 values in the 
interval curI £ (305, 310]. The dashed lines are interpolations 
between data points and merely serve as a guide to the eye. 
The other parameters are Uo = — IOojr and k = u;r. 


Fig. 8, which represents the essence of the underlying 
physics. For large negative cavity detunings Ac the 
photon-number distribution follows the mean values ob¬ 
tained for the lowest-band approximation from our semi- 
classical model [Eq. (14)[. With increasing Ac higher 
bands to -I- 2 appear when heating depletes the lowest 


































band if Eq. (18) is satisfied. Note that the population 
of higher momenta at relatively small photon number in 
Fig. 8 (e.g., at Ac = —T.Swr) can be traced back to 
trajectories like the one shown in Fig. 6d. 


V. CONCLUSIONS AND OUTLOOK 

We have shown that the dynamics of a quantum par¬ 
ticle trapped in a cavity-sustained optical lattice can 
reach a multitude of quasi-stationary solutions at the 
same operation parameters. Fluctuations of the cavity 
field as well as quantum dynamics of the particle even¬ 
tually trigger transitions between such states observable 
in single Monte Carlo quantum trajectories and stable 
for extended periods. Key properties of these strongly 
correlated atom-field solutions can be understood from 
an analysis in terms of localized Wannier functions and 
a mean-field approximation of the cavity mode. Quan¬ 
tum simulations exhibit few but fast transitions between 
such quasi-stationary states. Averaged over a sufficiently 
large ensemble, the final density matrix in this regime 
is a mixture of several Bloch bands with corresponding 
photon-number distributions. While the density matrix 
is mostly a mixture of such quasistationary states, some 
atom-field entanglement can be present in the transition 
phase, where the photon-number expectation value lies 
in between two stationary values. 
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Appendix A: Linear stability analysis of 
self-consistent solutions 


amplitudes via the bunching parameter b [Eq. (5)]. This 
dependence has to be considered in our stability analysis. 
We assume that 


a = ass + S, 


where 


and obtain 


d = 


S' 

S* ' ’ 


(A4) 

(AS) 


^ (q^ss + <5) — U{ass + S, a*g -|- (5*) (ass + S) + rj 


= C/(ass, + 


' dU[a, a*) 

■ 5 

dot 

OL — OLss 


ass + 0(<5'), 


(A6) 


where the term in the square parenthesis has to be inter¬ 
preted componentwise for each element of the matrix U. 
Using daf = {dnf) {°‘a) = a *we find 


' dU[a, a*) 

■ S 

da 

CX . — CX.SS 


as 


. 9Aeff 
* dn 


(a,V + ass5*) (J 


CXa 


= I 


dAeff 


dn 


Finally, we arrive at 


^^S = AS + OiS^), 


with the coefficient matrix 

.dAeff 


A := C/(ass, aL) +A 


dn 


d. (A7) 


(AS) 


■ (A9) 


We here present the derivation of Eq. (16). The equa¬ 
tions of motion for the field amplitudes a = {a,a*)'^ 
read 


Linear stability requires negative real parts of the 
eigenvalues of A. From 

Tr(A) = -2k < 0 (AlO) 


— a = U{a, a*)a + rj, 


(Al) 


with T] = (?7,??)^ and 

fiAes{a, a*) - K 0 \ 

^(“•“>=( 0 -iA.,(a, «•)-«)■ 


Its steady-state solution reads 


/-y - [ *^eff ^ 

'-^SS - \ — Tj 


(A2) 

(A3) 


We now investigate the stability of this steady-state so¬ 
lution against small perturbations. The effective detun¬ 
ing Aeff := Ac — Uob [cf. Eq. (17)] depends on the field 


follows, that the eigenvalues of A must be of the form 
-^ 1,2 = 0 ‘i ,2 i ib. At this point we have to discriminate 
between three cases: (i) Ai _2 = ai ,2 i ib, (ii) complex- 
conjugate eigenvalues Ai _2 = a i ib, and (iii) real eigen¬ 
values Ai ^2 = 01 , 2 - The first case can be ruled out since 
det(A) is real (see below). In the second case it follows 
from the negative trace [Eq. (A10)[ that a < 0, since 
Tr(A) = 2a, i.e., linear stability is always ensured. The 
third case implies det(A) = 0102 , which is positive for 

! 

negative eigenvalues. We thus require det(A) > 0, which 
evaluates to 


-k 


A 1 n 


1 '^SS 0 

on 

n—riss 


3Aeff 

dn 


> 0 . 

(All) 
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Inserting the steady-state photon number [cf. Eq. (A3)] 




(A12) 


simple arithmetic yields the condition 

dAcs 


K + Agfj -|- 2Aeff- 


K' 


_ 

2 + dn 


> 0, (A13) 


which can be recast into the form 


1-A 


dn -I- Agjj 


> 0 . 


(A14) 


This is precisely the inequality (16). 
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